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ABSTRACT 

We study tight bounds and fast algorithms for LCLMs of 
several linear differential operators with polynomial coef- 
ficients. We analyse the arithmetic complexity of existing 
algorithms for LCLMs, as well as the size of their outputs. 
We propose a new algorithm that recasts the LCLM compu- 
tation in a linear algebra problem on a polynomial matrix. 
This algorithm yields sharp bounds on the coefficient degrees 
of the LCLM, improving by one order of magnitude the best 
bounds obtained using previous algorithms. The complexity 
of the new algorithm is almost optimal, in the sense that it 
nearly matches the arithmetic size of the output. 

Categories and Subject Descriptors: 
1.1.2 [Computing Methodologies]: Symbolic and Alge- 
braic Manipulations — Algebraic Algorithms 

General Terms: Algorithms, Theory. 

Keywords: Algorithms, complexity, linear differential op- 
erators, common left multiples. 

1. INTRODUCTION 

The complexity of operations in the polynomial ring K[a;] 
over a field K has been intensively studied in the computer 
algebra literature. It is well established that polynomial 
multiplication is a commutative complexity yardstick, in the 
sense that the complexity of operations in K[a;] can be ex- 
pressed in terms of the cost of multiplication, and for most 
of them, in a quasi-linear way. 

Linear differential operators in the derivation d = and 
with coefficients in K(x) form a non-commutative ring, de- 
noted WL(x)(d), that shares many algebraic properties with 
the commutative ring K[x]. The structural analogy between 
polynomials and linear differential equations was discovered 
long ago by Libri and Brassinne [26, 8, 14]. They introduced 
the bases of a non-commutative elimination theory, by defin- 
ing the notions of greatest common right divisor (GCRD) 
and least common left multiple (LCLM) for differential op- 
erators, and by designing a Euclidean-type algorithm for 
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GCRDs. This was formalised by Ore [27, 28], who set up a 
common algebraic framework for polynomials and linear dif- 
ferential operators (and other skew polynomials, including 
difference and q-difference operators). Yet, the algorithmic 
study of linear differential operators is currently much less 
advanced than in the polynomial case. The cost of product 
in ¥L(x){d) has been addressed only recently in [22, 6]. 

The general aim of this work is to take a step towards a 
systematic study of the complexity of operations in ¥L(x) (d) . 
We promote the idea that (polynomial) matrix multiplica- 
tion may well become the common yardstick for measur- 
ing complexities in this non-commutative setting. The spe- 
cific goal of the present article is to illustrate this idea for 
LCLMs. We focus on LCLMs since several higher level al- 
gorithms rely crucially on the efficiency of this basic compu- 
tational primitive. For instance, algorithms for manipulat- 
ing D-finite functions represented by annihilating equations 
use common left multiples for performing addition [34, 32]. 
LCLMs of several operators are also needed as a basic task 
in various other higher-level algorithms [2, 24, 12]. Our ap- 
proach is based on using complexity analysis as a tool for 
algorithmic design, and on producing tight size bounds on 
the various objects involved in the algorithms. 

It is folklore that Ore's non-commutative Euclidean al- 
gorithm is computationally expensive; various other algo- 
rithms for computing common left multiples of two opera- 
tors have been proposed [21, 34, 32, 25, 5, 1, 23]. As opposed 
to Ore's algorithm, these alternatives have the common fea- 
ture that they reduce the problem of computing LCLMs to 
linear algebra. However, few complexity analyses [17, 18, 5, 
23] and performance comparisons [25, 1] are available. 

Main contributions. As a first contribution, we design a 
new algorithm for computing the LCLM of several opera- 
tors. It reduces the LCLM computation to a linear algebra 
problem on a polynomial matrix. The new algorithm can be 
viewed as an adaptation of Heffter's algorithm [21] to several 
operators. At the same time, we use modern linear-algebra 
algorithms [35, 37] to achieve a low arithmetic complexity. 
Our algorithm is similar in spirit to Grigoriev's method [20, 
§5] for computing GCRDs of several operators. 

Before stating more precisely our main results, we need 
the following conventions and notations, that will be used 
throughout the paper. All algorithms take as input lin- 
ear differential operators with polynomial coefficients, that 
is, belonging to WL[x](d), instead of rational function coeffi- 
cients. From the computational viewpoint, this is not too se- 
vere a restriction, since the rational coefficients case is easily 
reduced to that of polynomial coefficients, by normalisation. 
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Figure 1: Costs of various algorithms for the LCLM computation of k operators of bidegrees (d, r) in (x,d). 
Algorithms marked by a star (*) also compute cofactors for the same complexity. 



For Li, . . . , L k in K[x](d), we write LCLM(Li, . . . , Lk) for 
the primitive LCLM of Li, . . . , Lk in K[a;]{<9). We say that 
an operator L G K[x] (d) has bidegree at most (d, r) in (x, d) 
if it has order at most r and polynomial coefficients of degree 
at most d. The cost of our algorithms is measured by the 
number of arithmetic operations that they use in the base 
field K. The constant 6 £ [2, 3] stands for a feasible exponent 
for matrix multiplication over K (see definition in Section 2), 
and the soft-0 notation 0() indicates that polylogarithmic 
factors are neglected. Our main result is the following. 

Theorem 1 Let L\, . . . , Lk be operators in K[x](9) of bide- 
grees at most (d,r) in (x,d). Then LCLM(Li, ...,Lk) has 
order at most kr, degrees in x at most dk(rk — r + 1), and it 
can be computed in O(k 20 r 6 d) arithmetic operations in K. 

The upper bound dk(rk — r + 1) on coefficient degrees is 
sharp, in the sense that it is reached on generic inputs. It im- 
proves by one order of magnitude the best bound 0(k 2 r 2 d) 
obtained using previous algorithms. Moreover, for fixed k, 
the cost of the new algorithm is almost optimal, in the sense 
that it nearly matches the arithmetic size of the LCLM. 

As a second contribution, we analyse the worst-case arith- 
metic complexity of existing algorithms for LCLMs, as well 
as the size of their outputs. For instance, we show that the 
extension to several operators of the "folklore" algorithm 
in [34, 32] has complexity O(k B+1 r 0+1 d). We call this exten- 
sion van Hoeij's algorithm, after the name of the implemen- 
tor in Maple's package DEtools of one of its variants. These 
estimates are in accordance with our experiments showing 
that our new algorithm performs faster for large order r, 
while van Hoeij's algorithm is well suited for large k. 

Using our tight degree bounds, we also show that any 
algorithm that computes the LCLM of two operators of bi- 
degree (d, r) in (x, d) in complexity 0(r a d /} ) can be used as 
the building block of a divide-and-conquer (DAC) algorithm 
that computes the LCLM of k operators of bidegree (d, r) 
in complexity 0(k a+2l3 r a+ ^d 13 ). The costs of several algo- 
rithms are summarised in Figure 1, where notation A + 
DAC indicates that algorithm A is used in a DAC scheme. 

As a third contribution, we prove an upper bound B m 
2k(d + r) on the total degree in (x, d) of nonzero common 
left multiples (not necessarily of minimal order). This is a 
new instance of the philosophy, initiated in [7], of relaxing 
order minimality for linear differential operators, in order to 
achieve better arithmetic size. While, by Theorem 1, the 
total arithmetic size of the LCLM is typically k 3 r 2 d, there 
exist common left multiples of total size 4k 2 (d + r) 2 only. 

A fourth contribution is a fast Magma implementation 
that outperforms Magma's LCLM routine. Experimental 
results confirm that the practical complexity of the new al- 
gorithm behaves as predicted by our theoretical results. 

Last, but not least, we have undertaken an extensive bib- 
liographic search, which we now proceed to describe. 

Previous work. Libri [26] and Brassinne [8] (see also [14]) 
defined the notions of GCRD and LCLM of linear differ- 



ential operators, and sketched a Euclidean-type algorithm 
for GCRDs. Von Escherich [15] defined the related notion 
of differential resultant of two linear differential operators. 
Articles [8, 15] contain the embryo of an algorithm for the 
LCLM based on linear algebra; that algorithm was explic- 
itly stated by Heffter [21], and later rediscovered by Poole 
in his classical book [31]. The roots of a subresultant the- 
ory for differential operators are in Pierce's articles [29, 30]. 
Blumbcrg [3] gave one of the first systematic accounts of the 
algebraic properties of linear differential operators. Building 
on predecessors' works, Ore [27, 28] extended the Euclidean- 
type theory to the more general framework of skew polyno- 
mials. He showed [28, Theorem 8, §3] that, while the LCLM 
is not related to the GCRD by a simple formula as in the 
commutative case, there nevertheless exists a formula ex- 
pressing the LCLM in terms of the successive remainders 
in the Euclidean scheme. Almost simultaneously, Wedder- 
burn [40, §7-8] showed that the LCLM can also be computed 
by an extended version of the Euclidean-Ore algorithm, that 
computes Bezout cofactors along the way. 

In the computer algebra literature, algorithmic issues for 
skew polynomials emerged in the 1990s, and were popu- 
larised by Bronstein and Petkovsek [9, 10]. Grigoriev [20, 
§6] designed a fast algorithm for computing the GCRD of 
a family of linear differential operators; to do so, he proved 
tight bounds on the degree of the GCRD, by extending von 
Escherich's construction of the Sylvester matrix for two dif- 
ferential operators to an arbitrary number of operators. The 
bound is linear in the number of operators, in their maxi- 
mal order and in their maximal degree [20, Lemma 5.1]. 
Giesbrecht analysed the complexity of the LCLM compu- 
tation for two operators, but only in terms of their or- 
der [17, 18]. (Strictly speaking, his method was proposed 
for a different Ore ring, but it extends to more general set- 
tings, including the differential case.) For two operators 
L\,L2 € ¥L(x){d) of orders at most r, the first (Heffter-style) 
algorithm [17, Lemma 5] computes LCLM(Li,L2) in 0(r B ) 
operations in K(x), while the second one [18, Lemma 2.1] 
(based on the extended Euclidean-Ore scheme) uses 0(r 2 ) 
operations in K(x). To our knowledge, no algorithm cur- 
rently exists similar to the Lehmer-Knuth-Schonhage half- 
gcd algorithm [16, Chapter 11], using a number of oper- 
ations in K(x) that is quasi-linear in r. Li [25] pointed 
out that algorithms for the LCLM computation that have 
good complexity with respect to the order, such as the naive 
Euclidean-Ore algorithm, do not necessarily behave well be- 
cause of coefficient growth. He developed a generalisation 
of the classical subresultant theory to Ore polynomials, that 
provides determinantal formulas and degree bounds for the 
GCRD and the LCLM [25]. He also compared the practical 
efficiency of Maple implementations of several algorithms. 

Giesbrecht and Zhang [19, Theorem 2.1] mention a com- 
plexity bound of 0(r 5 d 2 ) for the LCLM computation of two 
operators of bidegree (d,r) in (x,d), based on an unpub- 
lished 2002 note of Li. Over fields of characteristic zero, 
Bostan [5, Chapter 10] sketched a general strategy for com- 



puting several constructions on differential operators (in- 
cluding LCLMs), based on an evaluation-interpolation ap- 
proach on power series solutions. He stated, without proofs, 
several degree bounds and complexity results. For two oper- 
ators L\,L2 of bidegree (n,n) in (x, d), he announced that 
using fast Hermite-Pade approximation for the interpola- 
tion step yields an algorithm that computes LCLM(Li, L2) 



0(r 



operations. The approach was enhanced by van 



der Hoeven [23], who showed that the costs of the basic op- 
erations on differential operators can be expressed in terms 
of the cost of multiplication in K[a;]{<9}, and proved the com- 
plexity bound 0(n e+2 ) stated without proof in [5, §10.5]. 

2. PRELIMINARIES 

Let (K, 8) be a differential field, that is, a field K equipped 
with an additive map 8 : K — !> K satisfying the Leibniz rule 
8(xy) = x8(y) + 5(x)y for all x, y € K. We denote by K [d; 8] 
the ring of linear differential operators over the differential 
field (K, 8). A nonzero element L in K[d; 8] is of the form 



L = a r d T + a r _i9' 1 + ■ 



+ 00, 



where a r , a, — 1, . . . , ao £ K with a r 7^ 0. We call r the 
order of L, and denote it by ord(L). The noncommutative 
ring K [d; 8] is a left (and right) principal ideal domain, for 
which a Euclidean algorithm exists [27, 28]. 

Let L, Li, . . . , Lk be nonzero elements in K[d; 8]. Then L 
is called a common left multiple (CLM) of Li, . . . ,Lk if L = 
QiLi = ■ ■ ■ = Q k L k for some Q±, . . . , Q k G K [d; 8}. A com- 
mon left multiple of the least order is called a least common 
left multiple (LCLM). Two LCLMs of L\, L k are K- 
linearly dependent. 

Our main focus is on the particular case K = K(:r), the 
field of rational functions with coefficients in K, and 8 = -jj^ , 
the usual derivation with respect to x. In this case, we use 
the notation K{x)(d) for K[d\ 8], and LCLM(Li, . . . , L k ) for 
the primitive LCLM of Li,...,L k in K[x](9), that is the 
LCLM of L\, . . . , L k computed in ¥L(x){d) and normalised 
in K[x] (d) with trivial content. However, in order to keep the 
mathematical exposition as independent as possible of any 
particular case, we stick to the more general setting K[d; 8] 
whenever we discuss mathematical properties and bird's-eye 
view descriptions of algorithms. 

Polynomial and matrix arithmetic. The cost of our al- 
gorithms will be measured by the number of field operations 
in K they use. To simplify the presentation, we assume that 
polynomials in K[a;]< n (1. e., of degree less than n in x) can be 
multiplied within 0(nlog(n) loglog(n)) = 0(n) operations 
in K, using the FFT-based algorithms in [33, 11]. Most ba- 
sic polynomial operations in K[:r]< n (division, extended gcd, 
interpolation, etc.) have cost 0(n) [16]. We suppose that 9 
is a feasible exponent for matrix multiplication over K, that 
is, a real constant 2 < 9 < 3, such that two n x n matrices 
with coefficients in K can be multiplied in time 0(n 9 ). The 
current tightest upper bound is 9 < 2.3727 [39], following 
work of Coppersmith and Winograd [13], and Stothers [38]. 

The following result, due to Storjohann and Villard [35, 
37], will be helpful to estimate complexities for solving lin- 
ear systems arising from LCLM computations. Note that 
this is currently the best complexity result on polynomial 
linear algebra. The probabilistic aspects of the algorithms 
described in this article are entirely inherited from it. 



Theorem 2 [35, 37] Let M be an m x n matrix with en- 
tries in 'K[x\ < d- The rank p of M can be computed together 
with m — p linearly independent polynomial elements in the 
left kernel of M within 0(mn p e ~ 2 d) operations in K by a 
(certified) randomised Las Vegas algorithm. 

Moreover, if m — n, then the determinant of M can be 
computed using 0(n° d) operations in K. 

3. LINEAR FORMULATION FOR 
COMMON LEFT MULTIPLES 

In order to connect the computation of common left mul- 
tiples with linear algebra, we introduce some more notation. 
For a nonnegative integer n, we denote by K [d; 8]< n the K- 
linear subspace of K [d; 8] consisting of all linear differential 
operators whose orders are at most n. Moreover, we define 
a if-linear bijection 



4> n ■ K[d; 8]< n 



K" 



Yli=o a '^ ^ (a n ,a n -i, 



■ ,a ) 



For a nonzero element P £ if[(9;<5] of order m, and for an 
integer n with n > m, we define the Sylvester-type matrix 



Sn(P) ~ 



( <t>n{8 
<t>n (0 



P) \ 

n. — 1 — m p\ 



\ MP) J 



The matrix S n (P) has n — m + 1 rows and n + 1 columns. 
In particular, S„(l) is the identity matrix of size n+ 1. This 
matrix enables one to express multiplication by P in K [d; 8] 
as a vector-matrix product. Precisely, for Q £ K [d; 8] <„_,„, 



MQP) = 0n-m(Q)S n {P). 



(1) 



Let Li, . . . , L k be nonzero elements in K[d; 8]. For n > 
maxKKt ord(Li), the matrix 



/ 5„(Li) 



M n := 



S n (L 2 ) 



\ 



S n (L k ) 



(2) 



V S„(-l) S n (-1) 
has (fc + l)(n+l)-Eti ord ( L 

rows and k(n+l) columns. 
The following theorem is the main result of this section. 

Theorem 3 Let L\, . . . , L k be elements in K[d; 8] \ {0} of 
orders r\ , . . . , r k , and let n > ord(LCLM(Li, . . . , L k )). 

(i) If L is a common left multiple of L\, . . . , L k such that 
ord(L) < n and L = Q\L\ = ■ • • = Q k L k , then the vec- 
tor (0 n _ ri (Qi), . . . ,(p n - rk (Q k ),(fi n (L)) belongs to the 
left kernel of the matrix M n defined in (2) . 

(it) If the vector (ui, . . . , Ufc, u) is a nonzero vector in the 
left kernel of M n , where Ui G K n+1 ~ Ti for i = l,...,k 
and u 6 K n+1 , then u is nonzero, and cpn 1 ^) « s a 
common left multiple of L\, . . . , L k with respective 
left cofactors (/)~^ ri (m), . . . , <t>~-r k ( u *0- 

(Hi) If p is the rank of M n , then 

ord (LCLM(Li, . . . , L k )) =p + £- = i r t - k(n + 1). 



Proof. Suppose that L = QiLi for 1 < i < k. By (1), 

4>n-ri(Qi)Sn(Li) = 4> n {L), 

which is equivalent to i^n-r-j (Qi)S n (Li) + <fi n {L)Sn{— 1) = 0. 
Therefore the vector ((f>„- ri (Qi), ■ • • , 4>n-r k (Qk), <t>n(L)) be- 
longs to the left kernel of M„. The first assertion is proved. 

Conversely, suppose that (ui, . . . ,Ufc, u) is a nonzero vec- 
tor in the left kernel of M„. Then UiS„ (Li) + uS„(-l) = 
for all i with 1 < i < k. It follows from (1) that 

ct>~- ri {ui)Li = ^(u). 

Thus, u is nonzero, for otherwise, since cf> n is an isomor- 
phism, all the Ui would be equal to zero. The second asser- 
tion follows. 

To prove the last assertion, we set L = LCLM(Li, . . . , L k ) 
and £ = ord(L). Assume further that L = QiLi for all i 
with 1 < i < k. Then L, dL, . . . , d n ~ l L are common left 
multiples of L\, . . . , Lk of orders at most n, and such that 

&>L = (ffQi} Li for all 1 < i < k and < j < n - I. 

By the first assertion, for < j < n — £, the vector 

vj = (4> n - ri (VQi) , . . . , <t> n - Tk (9 J O fc ) , <f>n (Vl)) 

belongs to the left kernel of M n . These vectors are K- 
linearly independent because L, dL, . . . , d n ~ l L are. On 
the other hand, if (ui, . . . , Ufc, u) is a nonzero vector in the 
left kernel of M n , where m € K n ~ ri+1 , ...,u fe € K n ~ rh+1 , 
and u € K n+1 , then ^(u) 

is a common left multiple of L\ , 
. . . , Lfc with order no greater than n by the second asser- 
tion. Hence, 0^ x (u) is a if-linear combination of d n ~ e L, 
d n ~ e ~ 1 L, . . . , L, because it is a left multiple of L. Hence, 
there exist c n -t, c n _^_x, . . . , Co in K such that 

^ 1 (u) = c n ~id n ~ l L + c n ~t~ 1 d n ~ t ~ 1 L + ■ ■ ■ + cqL, 

which implies that the last n + 1 coordinates of the vec- 
tor (ui, . . . , Ufe, u) — X^j=o c j v j are au ec l ua l to zero. Since 
this vector belongs to the left kernel of M n , all its coordi- 
nates are zero by the second assertion. We conclude that 
{vo, . . . , v n -i} is a Jf-basis of the left kernel of M n , and 
thus n — I + 1 is its dimension. Then (iii) follows from the 
rank-nullity theorem, because M n has (fc+1) (n+1) — X).^ 
rows. □ 

Since the rank of M„ is at most k(n + 1), a direct con- 
sequence of Theorem 3 (iii) is the following classical result. 

Corollary 4 For Li, . . . , L k g K[d; 5] \ {0}, 

ord (LCLM(Li, . . . , L k )) < ord(Li) + ■ ■ ■ + ord(i fe )- 

4. COMPUTING LCLMS 

In this section, we review a few known methods for com- 
puting LCLMs and present a new one based on Theorem 3. 

4.1 Computing an LCLM of two operators 

Given two nonzero elements L\ and L2 of respective or- 
ders ri and r-2, we consider various methods for computing 
their LCLMs. The first methods compute left cofactor(s) of 
the given operator(s) first, and find an LCLM by multipli- 
cation in K [d; 5]. The last method is specific to K = K(a;). 



4.1.1 Heffter 's algorithm 

The first method can be traced back to Brassinne [8], von 
Escherich [15] and Heffter [21]. The sequence: 

d^Li,... Li, d ri L 2 ,... ,8L 2 , L 2 

has ri+r 2 +2 elements, each of which is of order at most ri + 
r 2 . Thus, these elements are K- linearly dependent. To 
compute LCLM(Li, L2), the strategy is to find the maxi- 
mal integer m and corresponding elements ai,o, . . . , ai ;r2 _ m , 
a-2,0, ■ ■ ■ , a2,r!-m € K with ai, r2 _ m a2, ri -m / such that 

ai >r2 -id r2 ~ l L± + a 2 ,r 1 -jd ri ~ j L 2 = 0. 

i — m j=m 

Set Ai = EIi m ai,r2-^ r2 " l andA2 = EjL m a 2 , ri - 3 d ri ~ 3 . 
Then A\L\ + A 2 L 2 = 0. Therefore, the product A1L1 is an 
LCLM of L\ and L 2 due to the maximality of m. 

This method can be reformulated using the notation in- 
troduced in Section 3. For a vector v/0 represented by 

(0, . . . ,0,w fe+ i, . . .), where v k+1 / 0, 

k 

in a finite-dimensional A"-vector space equipped with the 
standard basis (1, 0, . . . , 0), (0, 1, 0, . . . , 0), . . . , (0, . . . , 0, 1), 
we define A/"(v) to be k. For n > max(ri, r 2 ), define 

= ( s st4 ) ■ (3) 

Then, Heffter's method consists in computing a vector v 7^ 
in the left kernel of U ri +r 2 such that A/"(v) is maximal. 

The next lemma connects the order of L = LCLM(Li, L 2 ) 
with the rank of U ri +r 2 - It easily follows from the observa- 
tion that a maximal subset of Jf-linearly independent ele- 
ments in {d T2 Li, . . . , Li, d ri L 2 , . . . , L 2 } consists of d T2 Li, 
. . . , Li and d e ~ r2 ~ 1 L 2 , . . . ,L 2 , where £ = ord(L). 

Lemma 5 Let Li, L 2 be two nonzero elements in K[d; 8] of 
orders r\,r 2 . Then ord (LCLM(Li, L 2 )) = rank((/ r i+r 2 ) ~ 1- 

4.1.2 Euclidean algorithms 

The second family of methods is based on the Euclidean- 
Ore algorithm for differential operators [28]. 

Ore's algorithm. Assume that ri > r 2 . Setting Ri = L%, 
Lti = L 2 , one can compute the Euclidean (right) divisions 

Ri ~ Ri-2 — QiRi-l, 

for quotients Qi £ K[d; 8], and remainders Ri ^ satisfying 
ord(Ri) < ord(Ri-i) for i = 3, . . . , m, and R m +i = 0. Then, 
as in the commutative case, R m is shown to be the GCRD 
of Li and L 2 . Ore [28, §2] proved that the following product 

Rm-lRm Rm-2Rm-\ ' ' ' R3R4 1 R 2 Rg 1 Rl (4) 

is an LCLM of L\ and L 2 . (Here AB~ X denotes the exact 
left quotient of A and B, that is Q such that A — QB.) 

Extended Euclidean-Ore algorithm. Wedderburn [40, 
§7-8] observed (see also [9]) that the computation of (4) can 
be avoided, if one replaces the Euclidean algorithm by its 
extended version. Precisely, letting Ci = 1, Ci = 0, and 

Ci = d- 2 - Qid-i, for i = 3,...,m, 

the product (C m _i — Q m -iC m )Ri is an LCLM of L\ and L 2 . 



Li's determinantal expression. As in the commutative 
case, a more efficient version of the extended Euclidean-Ore 
algorithm is based on subresultants [25, §5]. To avoid tech- 
nicalities, we present an alternative, efficient, variant of the 
subresultant algorithm, based on a determinantal formula- 
tion [25, Proposition 6.1]. This method assumes that the or- 
der g of the GCRD of Li and Li is already known. Then, one 
constructs a square matrix C of size ri+r 2 — 2g+2 whose first 
r\ + T2 — 2g + 1 columns are the first n + r 2 — 2g + 1 columns 

of the matrix ( 1+r2_s ^ 1 | j ; an d whose last column is 
the transpose of the vector (d T2 ~ a , . . . , d, 1, 0, 0, . . . , 0). 

ri—g+l 

If det(£) is denoted U, then ULi is an LCLM of Li and L 2 . 

4. 1.3 Van der Hoeven 's algorithm 

The algorithm that we very briefly mention now is specific 
to the case K — ¥L(x), where the base field IK has character- 
istic zero. It works by evaluation-interpolation. The idea, 
originating from [5], is to perform operations on differen- 
tial operators by working on their fundamental systems of 
solutions. Due to space limitations, and in view of its com- 
plexity analysis, we do not give more details here, and refer 
the reader to the article [23]. 

4.2 Computing an LCLM of several operators 

Given several nonzero operators Li, L%, . . . , Lk £ K[d; 8] 
we describe various ways to compute LCLM(Li, . . . , Lk). 

4. 2. 1 Iterative LCLMs 

An obvious method is to compute an LCLM of k operators 
iteratively, that is, 

L = LCLM (Li, LCLM(L 2 , . . . , LCLM(L fc _i, L k )) . (5) 

A computationally more efficient (though mathematically 
equivalent) method is by a divide-and-conquer algorithm, 
based on the repeated use of the formula 

L = LCLM(LCLM(Li, . . . , L Lfc/2J ), 

LCLM(L Lfc/ aj+i,.-.,ifc))- 

Of course, the efficiency of an iterative algorithm depends 
on that of the algorithm used for the LCLM of two operators. 
This is quantified precisely in Section 5. 

4.2.2 Van Hoeij's algorithm 

Another algorithm for computing the LCLM of k lin- 
ear differential operators was implemented by van Hoeij as 
Maple's DEtools [LCLM] command; it seemingly was never 
published. For k — 2, the method is folklore; it is implicit, 
for instance, in the proof of [34, Theorem 2.3]. A variant of 
it is also implemented by the 'diff eq+dif f eq f command in 
Salvy and Zimmermann's gfun package for Maple [32]. 

Informally speaking, the method consists in considering 
a generic solution hj of Lj for 1 < j < k, then in find- 
ing the first linear dependency between the row vectors 
d 1 (hi , . . . , hk ) = (d 1 ■ hi, . . . ,d l ■ hk). In order to per- 
form actual computations, these vectors are represented by 
the canonical forms (rem(9\ Li), . . . , rem(9 I , Lk)) , for i = 
0, 1, . . . , where rem(A, B) denotes the remainder of the right 
Euclidean division of A by B. Let 

L — a n d n + an-id™^ 1 + ■ ■ ■ + a , 



where ao, «Xi, . . . , a n are undetermined coefficients in K. For 
all i with 1 < i < k, let Ri be the right remainder of L 
in the division by Li. Then L = Ri mod Li. Since L has 
generic coefficients, ord(i?;) is equal to ri — 1, e.g., by [25, 
Lemma 2.3]. Note that ao, 01 , . . . , o n depend linearly on the 
coefficients of the Ri's. There are s = ri+" -+rk coefficients 
in ili, ■ • • , iife- Equating Ri = 0, for i = 1, . . . , k, we obtain 
a linear system 

(a n , a„-i, . . . , a )H„ = (0, 0, . . . , 0), 

where H n is an (n + 1) x s matrix over K. Thus, comput- 
ing LCLM(Li, . . . , Lk) amounts to computing a nontrivial 
vector v in the left kernel of H s with A/"(v) being maximal. 
The rank of H 3 is equal to the order of LCLM(Li, . . . , Lk), 
e.g., by [1, Proposition 4.3]. Note that the original version of 
van Hoeij's algorithm does not make use of this last fact, and 
potentially needs to solve more linear systems, thus being 
less efficient when the LCLM is not of maximal order. 

4.2.3 The new algorithm 

As a straightforward consequence of Theorem 3, the 
LCLM(Li, . . . , Lk) can be computed by determining a non- 
trivial vector v in the left kernel of M s given in equation (2), 
with A/"(v) being maximal. This method computes not only 
the LCLM, but also its left cofactors Qi, Q2, ■ ■ ■ , Qk, while 
van Hoeij's algorithm does not compute any cofactor. 

5. ALGORITHMS, BOUNDS, COMPLEXITY 

In this section, we let K = K(x) be the field of rational 
functions with coefficients in a field K, and S = ^ be the 
usual derivation with respect to x. Recall that in this case we 
use the notation K(x) (d) for K [d; S], and LCLM(Li, . . . , Lk) 
for the primitive LCLM of Li, . . . , L k in K[s](9). 

All algorithms analysed below are specialisations of the 
algorithms reviewed in the previous section to K = K(a;). 
Moreover, we make the non-restrictive assumption that all 
algorithms take as input linear differential operators with 
polynomial coefficients, that is, belonging to K[x](d). 

The degree of a nonzero operator L g K[s](9), denoted 
deg x (L), is defined as the maximal degree of its coefficients. 
As in the case of usual commutative polynomials, 

degJAB) = deg I ( J 4)+deg^(B) for all A,B £ K[x](d) \ {0}. 

5.1 Tight degree bounds for the LCLM 

First, we give a sharp degree bound for LCLMs. As we 
show later, this bound improves upon the bound that can 
be derived from van Hoeij's algorithm. 

Theorem 6 Let Li ,. . .,Lk be operators in K[x] (d) \{0}. Let 
s — ord(Li) + ■ ■ ■ + ord(Lft), and d — max*^ deg x (Li). 
i/L = LCLM(Li,...,L fe ), then deg x (L) < d(k(s + 1) - s). 

PROOF. By Corollary 4, ord(L) < s. It follows from The- 
orem 3 and Cramer's rule that every nonzero coefficient 
of L is a quotient of two minors of M s . Note that every 
square submatrix of M s has size at most k(s + 1), since M 3 
has k(s + 1) + 1 rows and k(s + 1) columns. Thus, the 
degree of the determinant of such a submatrix is bounded 
by d(k(s + 1) — s), because every entry of M s is of degree at 
most d, and the last s + 1 rows of M s are free of x. □ 
As a consequence of Corollary 4 and Theorem 6, the first 
part of Theorem 1 is easily deduced. 
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Construct the LCLM from K,. 
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Return the LCLM. 



Figure 2: Pseudo-code for Heffter's algorithm, van Hoeij's algorithm and our new algorithm. 



5.2 LCLMs of two operators 

The following result encapsulates complexity analyses of 
LCLM algorithms for two operators. Heffter's, van Hoeij's 
and our new algorithm are summarised in Figure 2. 

Theorem 7 Let L\,Li G K[x](9) be operators of bidegrees 
at most (d,r) in (x,d). Then it is possible to compute the 
LCLM of Li and L2 in complexity 

(a) 0(min(r e d 2 , r 3 d)) by Heffter's and van der Hoeven's 
algorithms, 

(b) 0(r B+1 d) by Li's and by van Hoeij's algorithms, 

(c) O(r d) by the new algorithm. 

PROOF. By [23, Theorems 5, 8 & 23], and using bounds from 
Theorem 6, the complexity of van der Hoeven's algorithm is 
0(min((rd) 2 r B ~ 2 , r 3 d)). The most costly parts of Heffter's 
algorithm are Steps 2, 4 and 6. Since the matrix U ri +r 2 has 
size 0(r) and polynomial coefficients of degree at most d, 
the rank and kernel computations involved in Steps 2 and 4 
can be performed using 0(r e d) operations, by Theorem 2. 
Step 6 consists in multiplying two operators in K[s]{<?) of 
bidegrees at most ((r — l)d,r) and (d, r) in (x,d). This can 
be done using 0(min((rd) 2 r 6 ~ 2 , r 3 d)). This proves (a). 

The dominant parts of Li's algorithm are the computa- 
tion of g = ord(GCRD(Li, L2)), and the expansion of O(r) 
minors of a polynomial matrix of size 0(r) and degree at 
most d. By using [20, Lemma 5.1] and Theorem 2, g can 
be computed using 0(r e d) operations in K, and the minors 
can be expanded in 0(r e+1 d). The dominant parts of van 
Hoeij's algorithm are Steps 2 and 4. Since k = 2, matrix H s 
has size 0(r). By an easy induction, its (r + j)th row has 
polynomial coefficients of degrees at most 2jd, thus H 3 has 
degree 0(dr). By Theorem 2, the rank and kernel compu- 
tations have complexity 0(r +1 d). This proves (b). 

The dominant parts of the new algorithm are Steps 2 
and 4. Since k = 2, the polynomial matrix M„ has size 
0(r) and degree at most d. By Theorem 2 again, the rank 
and kernel computations have cost 0(r B d). This completes 
the proof. □ 

Quite surprisingly, the costs of Heffter's and of van der 
Hoeven's algorithms are penalised by the complexity of mul- 
tiplication of operators, which is not well-understood yet for 



general bidegrees. Precisely, it is an open problem whether 
two operators of bidegree (d, r) in (x, 9) can be multiplied 
in nearly optimal time 0(r e ~~ 1 d). If such an algorithm were 
discovered, then the costs of both algorithms would become 
O(r d), improving the corresponding entries in Figure 1. 

5.3 LCLMs of several operators 

We analyse three algorithms for LCLMs of several opera- 
tors: DAC, van Hoeij's and our new algorithm. 

5.3.1 LCLMs by divide-and-conquer 

Theorem 8 Suppose that we are given an algorithm com- 
puting the LCLM of two differential operators which, on in- 
put Li,Z/2 G K[x](9) of bidegree at most (D,R) in (x,d), 
computes LCLM(Li,Z/2) in complexity 0(R a D^) for some 
constants a > 2 and f3 > 1 independent of D and R. 

There exists an algorithm which, on input Li, . . . , G 
K[x]{<9) of bidegrees at most (d,r) in (x,d), computes L = 
LCLM(Li, . ..,L k ) using 0(fc Q,+2/ V a+,3 d' 3 ) operations in K. 

Proof. Suppose without loss of generality that k — 2 e is a 
power of 2. To compute L, we use a strategy based on (6), 
similar to that of the subproduct tree [16, §10.1]: we parti- 
tion the family (Li, . . . , Lk) into pairs, compute the LCLM 
of each pair using algorithm A available for two operators, 
remove the polynomial content, then compute LCLMs of 
pairs, and so on. Let L[ a:b ] denote the LCLM of L a , . . . , Lt 
with the content removed. At level 1, the algorithm com- 
putes the k/2 operators Ln : 2], . . . ,L\k-i:k]: a ^ level 2 the 
fc/4 operators L[i : 4], . . . , L[ k _ 3:k ], and so on, the last compu- 
tation at level £ being that of L as the LCLM of L[i : fc/2] and 
L[k/2+i-.k]- Let C(R,D) — 0(R a D 13 ) denote the complexity 
of algorithm A on inputs of bidegrees at most (D,R). By 
Theorem 6, the operators computed at level 1 < j < £ have 
bidegree at most (2 j r,2 j d((2 j - l)r + 1)). Thus, the total 
cost of the DAC algorithm on k inputs of bidegree at most 
(d, r) is bounded by Y%Zl z3Tr ' C (2 J r, 2 3 d{(2 J - l)r + 1)) , 
plus the cost of the content removal, which is negligible. 
Up to polylogarithmic factors, the cost is bounded by 

E ^TT ' ( 2 ^r • (Vdrf = \ ■ r^ ■ d? ■ g (2T +^-i, 

3=0 3=0 

which is 0(fc a+2 'V +,3 d /3 ). This concludes the proof. □ 



The cost of the algorithm is essentially that of its last step; 
this is a typical feature of DAC algorithms. A similar anal- 
ysis shows that the iterative algorithm based on formula (5) 



is less efficient, and has complexity 0(k 



a + 2/3 + 1 a 



As a corollary of Theorems 7 and 8, we get a proof of the 
complexity estimates in the first three entries of Figure 1. 

5.3.2 Van Hoeij's and the new algorithm 

Theorem 9 Let L\, . . . ,Lk £ have bidegrees at 

most (d,r) in (x,d). One can compute LCLM(Li, . . . , Lk) 

(a) in 0(k e+1 r e+1 d) operations by van Hoeij's algorithm, 

(b) in 0(k 2e r 6 d) operations by the new algorithm. 

Proof. The proof is similar to that of Theorem 7(b) and (c). 
The most costly parts of van Hoeij's algorithm are Steps 2 
and 4. Matrix H 3 has size O(kr) and polynomial coefficients 
of degree 0(krd). By Theorem 2, the rank and kernel com- 



putations have complexity 0((kr) krd) = 0(k 



d). 



This proves (a). The dominant parts of the new algorithm 
are Steps 2 and 4. The polynomial matrix M s has size 
0(k 2 r) and degree at most d. By Theorem 2, the rank and 
kernel computations have cost 0((k 2 r) e d) — O(k 20 r d). □ 
As a corollary of Theorem 9, we get a proof of the com- 
plexity estimates in the last two entries of Figure 1. Note 
that Cramer's rule applied to the matrix H s analysed in the 
previous proof yields the bound 0(k 2 r 2 d) on the coefficient 
degrees of the LCLM. This bound is improved by Theorem f . 

6. SMALLER COMMON LEFT MULTIPLES 

Our approach to computing more common left multiples 
(CLMs), that are generally not of minimal order, but smaller 
in total arithmetic size than the LCLM, is similar to the 
linear-algebraic approach used in Section 3. However, in- 
stead of considering a matrix encoding the d l Lj, with poly- 
nomial coefficients, we turn our attention to a matrix en- 
coding the x ll d l2 Lj, with constant coefficients. 

Existence of smaller CLMs. The new building block to 
consider is, for an operator P in WL[x]{d) of total degree A 
in x and d, and an integer N > A, the ( N ~£ +2 ) x ( N + 2 ) 
matrix C'n(P) with scalar coefficients whose rows represent 
the operators of the form x %1 d l2 P for < i\ + 12 < N, 
in any fixed order, and whose columns are indexed by the 
monomials of total degree at most N, in any fixed order. 

Let L\, Lk be elements of K[x]{(?), with respective 
total degrees Ai, . . . , A^. For N > max{Ai, . . . , A^}, let 
M' N (Li, . . . , Lk) be the matrix 

/C N (Li) \ 



M' N = M' N (L ly ...,L k ) = 

V-/ ( « 2+2) 

This matrix has m(N) rows and n(N) columns, where 

iV + 2 



n(N) = k 



Assuming all Aj equal to a same value A, the matrix M' N 
certainly has a nontrivial left kernel when m(N) > n(N), 
) > (fc — 1) ( JV + 2 ) , which happens when 



that is when fc^-^ 2 



N > B for B = 



kA + 



^/4k(k - 1)A 2 + f - 3 



< 2fcA, 



where the approximation holds for large values of k or A. 
Using A < d + r yields the main result of this section. 

Theorem 10 Let Li, . . . , Lk be elements in ¥L[x] (d) \ {0} of 
orders at most r, and with coefficients of degrees at most d. 
There exist nonzero common left multiples of total degree < 
2k(d + r) in (x,d), and total arithmetic size 0(k 2 (d + r) 2 ). 

Algorithms for CLMs. A simple algorithm for comput- 
ing common left multiples of total degree B « 2k(d + r) 
in (x, d) is based on the left kernel computation of the 
scalar matrix M' B (L\, . . . ,Lu). This matrix has sizes of or- 
der kB 2 /2 « 2k z (d + r) 2 . The cost of the procedure is 



O(k i0 (d- 



it is dominated by the kernel computation. 



To simplify the discussion, we assume in the remain- 
ing of the section that L%, . . . ,Lk have bidegrees at most 
(d, r) — (n, n). Then, Theorem 10 implies that, while the 
LCLM has order at most kn and degrees at most k 2 n 2 , there 
exist common left multiples of order and degree at most 4fcn. 
However, computing such a small multiple by the previous 
algorithm of complexity O(k 3e n 20 ) is more costly than com- 
puting the LCLM by the last two algorithms in Figure 1. 

Here we briefly sketch a faster algorithm for computing a 
common left multiple of order and degree at most 4fcn, based 
on Hermite-Pade approximation [5, Chapter 10]. One de- 
termines series solutions of the Li at order 0(k 2 n 2 ), takes a 
random linear combination / of them, computes its first 4kn 
derivatives, and outputs a Hermite-Pade approximant of 
/,/',..., / (4fcn) of type (Akn, Akn). The dominant com- 



plexity is that of the Hermite-Pade step, 0(k e+1 i 



[36]. 



A Fast LCLM Heuristic. As an interesting consequence 
of this fast CLM computation, we deduce a very efficient 
heuristic for LCLMs, asymptotically faster than all algo- 
rithms in Figure 1. It proceeds in 3 steps: (i) compute O(l) 
CLMs of order and degree at most Akn; (ii) take two random 
linear combinations with coefficients in K[x\\ {Hi) return 
their GCRD. The dominating steps are (i) and (Hi). By us- 
ing Hermite-Pade approximants for step (i) and Grigoriev's 
algorithm [20] combined with Theorem 2 for step (Hi), the 
total complexity is O(k e+1 n 0+1 ). This is nearly optimal, in 
view of the LCLM size k 3 n 3 . However, we are not yet able 
to turn this heuristic into a fully proved algorithm. 

7. EXPERIMENTS 

We implemented 1 two variants of our new algorithm in 
Magma V2.16-7 [4] and compared them with Magma's built- 
in LCLM routine (command LeastCommonLef tMultiple). 

Some experimental results are summarised in Table 1. We 
take as input k = 2 random operators in ¥ p [x]{d), each of 
bidegree (d, r) = (n,n) in (x, d), where p is a medium-sized 
prime and n is of the form [2 j/2 ] , for 2 < j < 11. Column 
New gives timings for the first variant of the new algorithm, 
that uses Magma's built-in polynomial linear algebra solver 
(the Kernel routine), while column New+S gives timings for 
the second variant, based on our own high-level implemen- 
tation of Storjohann's high-order lifting algorithm [35] . Col- 
umn (iV, D) displays the size N and the degree D of the poly- 
nomial matrix dealt with by algorithms New and New+S. 
The dominating part of these algorithms is the left kernel 
computation for a polynomial matrix of size (N + 1) x N and 

All computer calculations were performed on a Quad-Core 
Intel Xeon X5160 processor at 3GHz, with 8GB of RAM. 
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Table 1: Timings (in sec.) for LCLMs of k = 2 ran- 
dom operators in ¥ p [x](d) of bidegrees (n,n) in (x, d). 



degree D. The most time consuming part of New+S con- 
sists in 0(log N) polynomial matrix multiplications of size N 
and degree D. To facilitate comparisons, column MM(N, D) 
shows the total time taken by 10 products of random poly- 
nomial matrices of size N and degree D over F p . Finally, 
column output size displays the total arithmetic size of the 
computed LCLM, that is, its number of coefficients in F p . 

Several conclusions can be drawn from Table 1. First, 
Magma's LCLM tool exhibits an exponential arithmetic 
complexity behaviour (when passing from bidegree (n, n) to 
(n+1, n+1), timings are multiplied by a factor close to 1.5), 
but it is relatively efficient for small input sizes. Both vari- 
ants of the new algorithm are faster for n > 10, and New+S 
gains a factor 65 for n — 23, and almost 1300 for n = 46. 

Second, timings in column New exhibit a practical com- 
plexity proportional to n , which is inherited from Magma's 
linear algebra solver on polynomial matrices. In contrast, 
New+S has a practical complexity proportional to n 3 ' 5 (but 
with a higher proportionality factor). This good behaviour, 
closer to the theoretical complexity O(n 0+1 ) predicted by 
Theorem 1, is inherited from Magma's very efficient polyno- 
mial matrix multiplication, through Storjohann's algorithm. 

Finally, timings in column New+S grow nearly linearly in 
the corresponding output sizes given in the last column, and 
these sizes match exactly the sharp bounds in Theorem 1. 
This experimentally confirms that size bounds and worst- 
case complexity analyses predicted by our theoretical results 
are reached in generic cases. 
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